First-order phase transitions in non-equilibrium systems: New perspectives 
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First-order irreversible phase transitions (IPT's) between an active regime and an absorbing state 
are studied in two models by means of both simulations and mean-field stability analysis. Hysteresis 
around coexistence is the result of the interplay of the length of the interface, its curvature and a 
memory effect related to the phase which is being removed. A controversy on the occurrence of 
scale-invariance is clarified, conciliating the behavior of IPT's with its reversible counterpart. 
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The study of far from equilibrium systems continues 
to attract great attention and has become challenging 
subject of interest for many areas of research in physics, 
chemistry, ecology, catalysis, economy, social sciences, 
etc. An intriguing feature of those systems is the 

occurrence of irreversible phase transitions (IPT's) be- 
tween an active regime and an absorbing state where the 
system becomes trapped. After the work of Ziff et. al. 

our understanding of second-order IPT's has expe- 
rienced a rapid growth since they have unambiguously 
been placed either in the directed percolation or parity 
conserving universality classes Due to the lack of 

experimental feedback on second order IPT's, the huge 
activity in the field is mainly of academic-theoretical in- 
terest. In spite of the existence of experimental evidence 
on systems undergoing first order IPT's [§, our under- 
standing of this field is far from being satisfactory. There 
are important controversies and many aspects still re- 
main unexplored. In fact, the claim that power-law be- 
havior could also hold for first order IPT's is certainly 
a puzzle [Q. A similar controversy has recently raised 
in the field of reversible transitions ||. However, power 
law behavior can be identified as finite-size effects which 
vanishes in the thermodynamic limit ||. Furthermore, 
the existence of hysteresis, which is a signature of first 
order transitions in equilibrium systems, has so far, not 
been explored in detail in the field of IPT's. 

The aim of this work is to present an study of first 
order IPT's based on two complementary techniques, 
namely extensive numerical simulations and mean field 
approaches. Two different model systems, namely the 
ZGB model for a catalyzed reaction ||]4| and stochastic 
game of life (SGL) for a society of living individuals ||] 
are investigated in order to allow the discussion of useful 
comparisons. 

The SGL model is a cellular automaton defined on a 
square lattice where each site <7y can take only two val- 
ues <7y = {0 (dead site), 1 (living site)} and interacts 
with its eight nearest neighbors which defines its neigh- 
borhood (NH) . The system evolves from a given time to 
the next time step as follows: i) A living site whose NH is 
dead or allocates one living site, will die. ii) A living site 
whose NH allocates more than three living sites, will die. 
iii) A living site whose NH allocates two or three living 



sites, will survive with a probability p s . iv) A dead site 
whose NH is dead or allocates one living site, will remain 
in this state, v) A dead site whose NH allocates more 
than three living sites, will remain in this state, vi) A 
dead site whose NH allocates two living sites, will become 
a living site with a probability pb- vii) A dead site whose 
NH allocates three living sites, will become a living site. 

The phase diagram of the SGL has two phases, namely 
extinction and life, both separated by a first order coex- 
istence curve ||. A mean field (MF) equation for the 
time evolution of the density of living sites x is given by 
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where x is the order parameter and y = 1 — x. The 
fixed points of eq. (1) satisfy 4| = f(x) = 0, and the 
stable ones correspond to the stationary states of the 
system. We define a potential V(x) through the rela- 
tion = f{x) = -^M. Then, the stable (unsta- 



ble) fixed points correspond to the minima (maxima) of 
V(x), respectively. The values {pf ex ,pl° ex ) that satisfy 

^fci = and d2 y[ x) = determine the MF coexistence 

ax dx A 

curve. Figure 1 shows a plot of V(x) as a function of 
x and pb, keeping p s = 0.1 fixed, x = is always a 
minimum of V(x) which corresponds to the absorbing 
state. In addition, the surface exhibits another valley at 
a higher density which corresponds to the stationary liv- 
ing phase. The valley at a higher density is vanishing 
when decreasing pb, and it actually disappears at the co- 
existence point pc° ex given by (p s = Q.l,pl° ex = 0.3131) 
(see figure 1). Further decreasing pb leads to potential 
functions displaying only the valley at x = 0. Figure 1 
clearly shows that the valley at a higher density disap- 
pears at a value x c well above x — 0. So, a sharp jump 
in the order parameter of the system is observed indicat- 
ing a first order IPT in agreement with the simulation 
results gj. Since the system always evolves to the po- 
tential valleys, the stationary state will depend upon the 
initial density. This MF result is also in agreement with 
simulations. 

The ZGB model is an approach to the catalytic ox- 
idation of carbon monoxide, CO + l/20 2 -> C0 2 [@. 
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The arrival probabilities of CO and Oi are normalized, 
Pco + Po 2 = 1, so one has a single parameter which 
is selected to be Pco- The ZGB model exhibits a first 
order IPT between a regime of sustained reaction and a 
CO-poisoned state close to P 2 co = 0.5256(1) @. For 
details on the ZGB model see e.g. [pPJll[|. 



ter has occupied most of the system. Since the density 
can only grow from the interface between the coexisting 
phases by invading the absorbing cluster, a larger value 
of pb is needed. This is observed in figure 2, where the 
increasing density branch is the rightmost branch of the 
loops. 
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FIG. 1. 3d plot of the potential V(x) for p 3 = 0.10, corre- 
sponding to the SGL model. Darkest areas correspond to the 
potential minima. 



In this paper, hysteresis in the SGL model has been 
studied using the constant-coverage (CC) ensemble @- 
First, one has to achieve a stationary configuration using 
the standard ensemble. Then, the system is switched to 
the CC ensemble where the density x of living sites is 
varied stepwise. Since the SGL model has two external 
parameters, p s and pb, in the CC ensemble p s is kept 
constant and the control parameter is now x. Each point 
of the phase diagram, x versus pb (see figure 2), is eval- 
uated after allowing the system to relax tr time steps. 
The phase diagram is identical to that obtained with the 
standard method only within the living phase, i. e. for 
Pb > Pb° ex ■ When approaching the coexistence point, the 
curve bends downwards, close to pl° ex , becoming par- 
allel to the x axis. Then, the value of the parameter 
Pb remains constant, within this density range, pointing 
out that the living and the absorbing phase coexist. The 
MF phase diagram (x versus pb) exhibits two branches of 
fixed points, the upper stable curve and the lower unsta- 
ble one, that intersect at the coexistence point where the 
sharp jump in the density occurs. Further decreasing the 
parameter pb, leads to only one stable stationary state, 
namely the absorbing state. Then, no stable stationary 
living state exists for densities < x < x c . Using the the 
CC ensemble the system is forced to support a density 
within this range, so it splits in two stable phases, namely 
the living phase at pb = pg oe:E and the absorbing phase 
that behave almost independently. This is the reason 
why pb remains constant in the coexisting region. After 
achieving a density close to the absorbing state, we begin 
to slowly increase it. At that stage, a big absorbing clus- 
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FIG. 2. Hysteresis loops obtained by means of the CC en- 
semble for the SGL model, keeping the parameter p s =0.1 
fixed, for different relaxation times tr. Loops are generated 
counterclockwise (see arrows). The dashed line indicate the 
position of the coexistence point. 

This invading process, which can be regarded as a 
memory effect related to the absorbing cluster, is unique 
for the growing branch. It is also found that the longer 
the relaxation time 77? the narrower the loops. In addi- 
tion, the left branch (decreasing density branch) seems 
to fall downwards at the same pb value independent of 
the value of tr. Then, the right branch approaches the 
left one for longer tr. After a view inspection of two con- 
figurations obtained at the same density but on different 
branches (not shown here), it is clear that the presence of 
a large absorbing cluster in the configuration correspond- 
ing to the right branch and the existence of many small 
clusters in the ones corresponding to the left branch is 
the main difference among them. Then, the length of the 
interface between the two phases is clearly larger for con- 
figurations on the left branch, and consequently they are 
more efficient at keeping the population constant. This 
results in a lower value of the parameter pb- Notice that 
along the decreasing branch the memory effect is due 
to the presence of small active clusters. Also, loops be- 
come narrower for longer tr because the system evolves 
towards the most disordered state with the largest inter- 
face at a given density. 

In the case of the ZGB model, the density of CO 
{&Co) plays the role of the control parameter (figure 3). 
Hysteresis effects are absent in small lattices (L < 64), 
while hysteresis loops can be distinguished for larger 
lattices (figure 3). So, for the CO-growing (decreas- 
ing) branches of the loops we obtain P^ = 0.52641(1) 
(Pq = 0.52467(3)), respectively. Notice that the value 
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of the parameter remains almost constant for a wide 
range of C o values (0.1 < 9 CO < 0.90) for L > 512. 
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FIG. 3. Plots of 9co versus Pco obtained using the CC 
ensemble with tr = 100 mcs and lattices of different size. 
The point S shows the position of the upper spinodal point 
Pco- Arrows pointing up and down show the growing and 
decreasing 8co branches of the hysteresis loop, respectively. 



Assuming that the coexistence point lies in the middle 
of the loop, we get P2C0 = 0.52554(4). We claim that 
this value, obtained by measuring 400 different points 
each of them averaged over 2a;10 3 Monte Carlo steps, is 
the most accurate available |lC| . Snapshot configurations 
were obtained for both branches (not shown here). For 
the growing branch a massive CO— cluster preveals and 
after becoming large enough, it percolates along only one 
direction of the lattice, forming two (relatively flat) in- 
terfaces at the boundary of the coexisting phases. When 
6co is increased the interface roughens and eventually 
dangling ends evolving from both parallel interfaces get 
in contact causing the onset of percolation in both direc- 
tions of the lattice. At this point the system jumps to the 
other branch of the hysteresis loop. So, the interface be- 
tween coexisting phases has become shorter and its cur- 
vature has decreased. Smaller local curvature leads to a 
less efficient O2 adsorption mechanism and consequently 
a larger 2 pressure (smaller Pco) is needed in order to 
decrease the 9co- In addition, the memory effect asso- 
ciated to a bulky active phase, results in a bigger Pco 
value along the growing branch, in contrast to the mem- 
ory effect along the decreasing branch which is related 
to a bulky CO cluster. The phase diagram of the ZGB 
model as obtained by the CC-ensemble, clearly exhibits 
the upper-spinodal point P$o^ which strongly depends 
on the lattice size. In the infinite-size limit (L = 00), 
we obtain the value P^, (L = 00) = 0.5270(5) which 
may be compared with Pco(L = ^56) ~ 0.527 [pTj and 



P$o(L 



=?) ~ 0.5285 0. An stability analysis, simi- 
lar to that described above for the SGL, has also been 
performed within the MF approach proposed by Dick- 
man g]. Both models display two different branches of 
fixed points, namely a stable branch and an unstable one, 



which coalesce at a point identified as the spinodal point 
H . CC simulations clearly show the existence of a spin- 
odal point for the ZGB model but no evidence is found 
in the case of the SGL model. 

Another approach for the study of IPT's is the epi- 
demic analysis (EA) 0. For the SGL, we have studied 
the time evolution of the average number of living sites 
(N(t)) as it is shown in figure 4. Initializing the simu- 
lation with a small colony of living sites in an otherwise 
extinct state, it is found that the asymptotic regime is 
reached after t > 10 4 updates and it is possible to iden- 
tify subcritical, critical and supercritical curves as well. 
The critical point (p s — 0.1, pb — 0.47188) is in excel- 
lent agreement with the coexistence point obtained by 
means of the CC ensemble (see figure 2). Notice the high 
sensitivity of N(i) to tiny changes in the values of the pa- 
rameter. N(t) also displays a short time regime (t < 10 2 ) 
and an intermediate time regime (10 2 < t < 10 4 ), the last 
resembling a plateau behavior. 
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FIG. 4. Log-log plots of the average number of living sites 
N(t) versus t for different values of the parameter pb, keep- 
ing p s = 0.1 fixed. The EA are initialized using small living 
colonies (bigger colonies) for the main plot (inset plot), re- 
spectively. 

Starting with bigger colonies (inset of figure 4), we ob- 
serve the following three main differences, i) the short 
time behavior is absent, ii) intermediate time regime 
dominates from the very beginning and, iii) asymptotic 
regime is achieved sooner (t « 10 3 ). These features can 
be understood on the basis of the potential V(x). In fact, 
Figure 1 shows that for low initial densities, x will flow 
to x(t — > 00) = 0. This is the reason why N(t) displays a 
decreasing short time behavior in figure 4, which appears 
to be a universal feature observed in every EA of a first- 
order IPT. However, in some few cases and due to a den- 
sity fluctuation, the system may overcome the potential 
barrier, reaching a small region where V(x) « constant. 
Consequently, the system can remain for a long time in 
a region where x « constant until another fluctuation 
would drive it either towards x = or to the station- 
ary living state. This explains the plateau observed in 
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the intermediate time regime and the asymptotic regime 
as well. Another way of corroborating the above expla- 
nation is to initialize the simulation with larger colonies 
(see inset of figure 4). So, as the density is in a region 
where V(x) w constant from the very beginning, the 
short time behavior is absent and only the plateau and 
the asymptotic time behavior remain. The short time be- 
havior of the system can be easily confused with a power- 
law behavior, mainly if small lattices and short times are 
considered. 

In the case of the ZGB model the absorbing state is a 
lattice fully covered by CO. So, the EA is started with a 
covered lattice except by a small amount of empty sites 
placed at the center of the sample. Figure 5 shows results 
for the number of empty sites (N(i)) versus the time t. 
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FIG. 5. Log-Log plots of the number of vacant sites 
N(t) versus t for EA of the ZGB model. Results aver- 
aged over 5x10 s different runs (▼ P§ , ■ P§ Q , □ P 2C o, 
V Pco = 0.52345). For the later, two straight lines have 
been drawn for the sake of comparison: the dashed one with 
^eff — 2 anc i the full one with slope 2, respectively. The inset 
shows a semi-logarithmic plot of N(t)(T/t)~ 2 versus t/T with 
T = 183, according to equation 2. 

Results obtained for Pco < P§o show pronounced 
curvature, with clear evidences of a cut-off. So, the dy- 
namical critical behavior of the ZGB model at coexis- 
tence does not exhibit scale-invariance. This finding is 
in contrast with previous results claiming power-law and 
scaling behavior 0. For Pco > ^co' ^(*) exhibits a 
pseudo-power-law behavior over many decades with an 
effective exponent 7y e // — 2.0 ± 0.1. However, eventu- 
ally after a long time, a successful epidemic spreading 
preveals and N(t) suddenly growths as N(t) oc t 2 , as it 
is shown in figure 5, indicating a spatially homogeneous 
spreading. In order to understand the critical behavior 
of N(t) at coexistence we propose the following anzats 

N(t) ~ (t/T)-*" exp(-t/T), (2) 

where T sets a characteristic time scale which is validated 
by the plot in the inset of figure 5. A regression analysis 



gives T = 183 ± 3. Summing up, one observes a pseudo 
power-law behavior for short time (t < T) which crosses 
over to an asymptotic exponential decay for larger time. 

In summary, we have studied two different models that 
exhibit first-order IPT's, using extensive computer sim- 
ulations and MF approximations. Despite of the lack of 
universality, the critical behavior of both models share 
several features, namely (i) hysteresis is absent for small 
lattices, (ii) hysteresis effects are due to the interplay 
between differences in length and curvature of the inter- 
faces among the coexisting phases at the same density, 
and the existence of memory effects related to the phase 
which is being removed. Which of the three effects play 
a major role may depend on the system and is the sub- 
ject of an ongoing investigation, and (iii) the occurrence 
of power-laws in the dynamical critical behavior of first- 
order IPT's can be safely ruled out. This last finding 
conciliates the behavior of first order IPT's with their 
counterpart in equilibrium systems where it is well es- 
tablished that the existence of short range correlations 
inhibits the observation of scale invariance. However, 
major differences to be noted are: (i) in contrast to the 
ZGB model, spinodal points are not observed in the SGL 
model, (ii) for tr — > oo, hysteresis loops collapse to a sin- 
gle curve. However, for the SGL model there is a single 
vertical line at coexistence, while for the ZGB model the 
curve splits out in two vertical lines resembling a phase 
transition in the parameter space, and (iii) the interme- 
diate time dynamic critical behavior is different. 
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